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ABSTRACT 

Handling nonlinear systems 15 not a simple task, 
especially if a computer is not available. Generally, 
although a precise result for every case would be preferable, 
an approximate solution is usually sufficient for nase 
engineering solutions. 

R. Boxer and S. Thaler developed a technique [1] 
presented in Proceedings of IRE that allow o obtain the 
response of linear and nonlinear systems, without knowledge 
of the roots of the system characteristic equation. 

The aim of this work is to test this technique in 
Ши де “Пе solution of several time varying linear 
equations and to produce a computer program tnat improves 


the solution accuracy and ease of use. 
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I. INTRODUCTION 


Two basic methods exist for solving the response of a 
system which generally is presented as an integro-differential 
equation: analytic methods and numerical methods. 

For simple problems, analytic methods give precise 
solutions and have the advantage that these solutions can 
be given as a function of the parameters of the system. 

As the order of the differential equation increases, 
numerical methods give less labor in providing a solution 
and can be, at some extent the only method available. 

Laplace transform is limited to the case of linear 
differential equations. Nonlinear differential equations 
are out of the scope of this technique, preventing its 
application. 

Even sometimes with linear systems, if the inverse 
integral is complex, it is difficult and sometimes impossible 
to find analytically the time function. 

However, it is possible to approximate the Laplace 
transform of a linear system by a z-transform, obtaining 
thereafter a time series representing the values of the 
response at equally spaced intervals of time, by synthetic 
division. 

The method [1] deseribed and used in this paper utilizes 


the z-transform directly in order to obtain the solution of 


une Пъдезстнана nonlinear systems. 





In this paper, an explanation is given of this method 
applied to linear systems and its extension to time-varying 
and nonlinear systens. 

A computer was used for faster and more reliable results, 
but the method described is compatible with handwork. 

О ЈЕО. time по the computer, Гог this technique, is 
much smaller than other numerical methods due to the less 
number of iteration steps needed. This technique could 
reduce solution times by several orders of magnitude over 


conventional numerical methods. 





eee hoc AL LON OF THE METHOD 


A. PROCEDURE 
The basic procedures of the method are given, before 
the description or the theory, to show its simplicity. 
Suppose that one has the Laplace transform of an output 


response. It will be usually in the following form: 


а +а 5 +а 5° + а s£ + .., а ar 
_ © 1 2 К n 
F(s) = >, (1) 
bg t b,8 05 To... + 0,5 + eee + 0,5 


The Steps to obtain the time response are: 


1, Express the function F(s) as a rational fraction in 


powers of st by dividing the numerator and denominator 


by за 


"Substitute for sf a rational fraction in powers of 271 


obtained from Table I and rearrange F(s) as a rational 
fraction in powers of zt, 

3. Divide the resulting expression by T, where T is the 
time interval between points at which the result is desired. 


4. Expand the fraction by synthetic division into a series 


of the form: 


zu “г =n 
р t D; z + D5z к... + D, 2 + eee (2) 


where р, the coefficient of 27, is the approximate value 


of the time response at t = nT. 





B. LINEAR SYSTEMS 
Take an arbitrary function f(t) as in Fig. l(a). Its 


Laplace transform and inverse Laplace transform are defined 


аз: 


со 


Fis) =f f(t) e >” at (3) 
0 
c+ je 
f(t) = 3H го ее е5© аз (n) 
єг 





; Aca: 400) 
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4 
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fH 

: ПАЛКЕ an 
DNUS Arbitrary T(t) Fig. 1(b): Sampled f(t) 


Figure 1(b) shows f(t) sampled by shifted unit impulses. 


The result is a series of impulses at t = nT. 





Applying Laplace to the sampled function one gets: 


БК б “з v (at) в ni 


n=0 


er 


Substituting in (5) сс by z so: 


Z = E | (6) 
the result obtained is known as the z-transform: 


F*(z) = В ео (7) 


l t4 8 
e 


[1 


and the inverse is given by: 


Мат) - 6 Е%(2) 2071 az (8) 


иене но convour 15 a circle centered at the origin and 
medange ali the singularities of F(z). 

Е uae onn O) maps the 10 axis of the s-plane into a 
unit circle around the origin in the z-plane. 

There are other methods for obtaining direct and inverse 
z-transforns. 


In closed form the function will appear as: 


D = 
о 2 


F*(z) - z (9) 


10 





Division of the numerator by the denominator, will 
result in the coefficients of f(nT) as defined in (7), 
without the use of the inversion integral formula (8). 

In order to apply this useful characteristic of 
z-transform theory to non-sampled systems one uses the 
inverse Laplace integral (4), 

Оте 


f(t) = i i F(s) e5? as 
2T j 2-19 





which can be expressed piecewise as: 


min J ° 
f(t) = I d^ ee T | : (ауе rr de 
2011 JT ZP" 


(10) 


The way the contour is divided is shown in Fig. 2. The 
first integral gives the part of the contour shown between 
-(n/T) < w< n/T . The remaining contour w > (т/Т),ш < -(m/T) 
is given by the second integral. 

If 1/T is sufficiently large, a good approximation will 
be held by the first integral of (10) and the second 


integral may be discarded, being the resulting error given 


Буг 
Jo 
e(t) = 513 f  f(s)e9" 4 Е(-ва)е79" ав (11) 
= J HZ 


I 








Figure 2: Contour in s-plane for approximate solution 


The approximation, f,\t), on ft) is given by: 


Jl ZP 


MES st 
f,(t) 273 ns F(s) e ds (12) 
Making t - nT 
UTE 
Ё.(пЕ) y e (16) 
А 2T j . 


Solving eq. (6) for s, one obtains: 


SN SEN (14) 


ша 





Making the substitution of (14) into (13) the following 


equation is obtained 


ЕІ In/T 
РАСТ) = эт [| — p(1/T 1n z) ePT(17D)1n 2 ati5m 1а г) 
— /T 
(15) 
г. (пт) + ДО в 43 Б(1/Т 10:52) 20-1 ағ (16) 
А aTj ЈЕ 


where the contour, again, is the unit circle in the z-plane. 

If eq. (16) is compared with eq. (8), it is seen that it 
has the form of an inverse z-transform except for the constant 
multiplier 1/T. 

On the other hand, as F(1/T ln z) is a transcendental 
function, it is not expandable directly into a series of 
powers of zt by synthetic division. The synthetic division 
technique can, however, be used if a proper approximation 
courd be found. 


If F(s) is expanded in a series of descending powers of 


s like 
2 2 = d. 
F(s) = ав + ans * * (211) 
a solution of the form 
+2 
e) = ay + a,t + > 57 t ... (18) 


eensbesobcrained, but this solution is not always valid. 


15 





The way to follow is to rearrange F (1/T)/ln z asa 
1! 


rational function of two polynomials in powers of s ~ = T/ln z. 


For this purpose, the power series for Іп z used is: 
In z = 2[u 4 1/3-43 4 1/5 u” +...) (19) 


where: 


= l] -= 271 
1 


1+ 2 


So the following equation is obtained: 


1 T T/2 
ci s mE _ (20) 
= In z u + 1/3 u3 + 1⁄5 y? 


By synthetic division the following Laurent series is 


obtained: 
1. 6 1 h 3 11 54... 
A MO er 


To obtain series expansion of ээ bothsides of (21) 
can be raised to the gen power. 

Table I (on page 15) was obtained by using only the 
principal part and the constant term of the Laurent series. 


The result is of the form 


14 





TABLE I 


-1, 


Е (2 
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- | _ 
-K Nk (Е )-p (а (22) 


a The right hand 


k 


where N, is a polynomial in powers of zu 


side of (24) is called the "z-form" for s 


Example 1 — Third Order System (Fig. 3) 


221 
се 






O;(): £ 





Fig. 3: Typical Feedback System for Example 1 


With a step input, as shown, the Laplace transform of 


the output 15: 


1 
086 (s) = —— — n (23) 
0 oe 


Following the step by step procedure as given in the 


Introduction, 


99s) = — — (24) 


16 





Перо сиса the correspondent forms of s É from 


Table I and dividing the result by T: 


a + + 272) 


0 {жй =ош — Е 2777777 тан 
ын (12+6Т+Т°) — (36*6T-01^)z l + (36-6T-91^)z ^ — (32-6T«1^)z 2 


25) 


The solution is obtained by choosing T, the time interval 
at which rate our solution is obtained, and proceeding in 
dividing the denominator into the numerator. 


Let us assume T = 1,0. resulting in the following output: 


-1 -2 
Өөд" (2) = = = = (26) 
19-332 — + 21z * - Tz 
By long division, the result is: 
1 


The accuracy improves if smaller value for T is chosen. 


Letting T = 0.5 in this same example, substituting on eq. (2506 


1.521 4 1,57"? (27) 


Ө *(z) = ЕЕ к ЕЕ 
OA 15.25 - 36.752 = + 30.7527 – 9.252 


а. out the long division, the result is: 


1 15) 22: 06. 06120102 


The points obtained are plotted together with the precise 


solution in Figs. 4. 


il 





915 
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о X 
2 
©, 
O 
A l  _ l llelllJl._CCeCCYC——— ——— n = 
E 285 2) 034 0:5 225 019 
Figure 4 
к==еаје = 2.0 units/inch : . 
Е : Е — Precise solution 
Y-scale = ‚5 units/inch Y 2 


x Approximate Solution 
ШОР 2027-2100 

+ Approximate solution 
(Ol AL 5 


Example 2 — Solution of Linear Differential Equation 


To use the present technique in obtaining the solution 
of integro-differential equations, one takes the Laplace 
Ви nro tae equal ton. initial conditions included as 
usual, and then proceeds as before. However, after trans- 


forming the equation, the following rules can be applied: 


18 





i. If the degree of the denominator of the Laplace 
transform is higher than the degree of the numerator by two 
premore, proceed directly, 

2. If the degree of the numerator is one less, equal to, 
or greater than the degree of the denominator, divide until 
the remainder meets the criterion of Rule 1. 

Using the procedure stated in the rules above one 
eliminates the high frequency from F(s), that is implied if 
the numerator is one degree less than the denominator. 


To illustrate this, a simple example is shown: 


ax " . 
тг - (28) 
Assuming for initial condition Y(0) = u, the Laplace 


transform Y(s) 15: 


+ в. -u + 
Д SU _ SU 1 (29) 


О 723: 


Using Rule 2 above, one gets: 


E l-u 
Ir re i 


In this procedure, the step funetion of amplitude u is 
added to the result obtained after operating upon 


(1-w/(s**s). 


15 





Now, letting u = 0.1, and expressing the second term 


ОҒ еа. (30) in powers of 571, 


-2 
5.25... (31) 
ТЕЗ 


Y (s) = 
Substituting the z-forms of Table I (on page 15) and 


dividing by T, 


1 2 


Mid ë 0.9T(1+10z +z ^) T 


(1246T) - 2421 + (12 - 6T)z ^ 


Letting IT = 0.1 and carrying out the eng division one 


Бере: 
-1 -2 -3 -l 
О СО 2.2337 P+ 2062 "sve (33) 


This result is plotted in Fig. 5 after adding the step 


function of amplitude 0.1. 





о 
a 
# 
Е) + 
э 2 322 025 038 о 00 
Figure 5 
X-scale = 2.0 units/inch — Precise solution 
Y-scale = .5 units/inch + Approximate solution 
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С. TIME-VARYING SYSTEMS 
A linear differential equation with time-varying 


coefficients may be expressed as 


ап 80-1 
ES БН ES na cox] и 22: гт) 


ас? dt 


where f£ (t), В (t),... f. (t) are functions of the 


n-1 
independent variable t. 

To solve this type of equation the following Steps may 
be executed: 

1. Obtain the Laplace transform of (34) in the usual 
manner including initial conditions. 

ee) Proceed as іп the constant coefficient case discussed 
above. TIn performing the long division process, however, 
change the values of Cas...» at each step in the division 


process tothe values at the corresponding times of' the 


functions 


Ca = 221) 
(35) 
Су = fg(t) 
In general 
C, = f(t) (36) 


21 





D. NONLINEAR SYSTEMS 

A similar manner like used in II.C. may be applied to 
nonlinear differential equations. The regression equation 
is obtained by long division process as in the time varying 
example. 

However, the factors Са Cs os Со must now 
represent functions of the dependent variable. 

Generally, a nonlinear differential equation can be 


represented as: 


n 1-1 
а а dY 
E eg peor] и О на 
at? i ach | 1-1 dt? 
+ EME EOS (37) 


Perne собе rules tor the procedure as follow: 

1. Obtain the Laplace transform of (37) in the usual 
Шээс ас initial conditions. 

2. Proceed as in the constant coefficient case. In 
performing the long division process, however, change the 
value of Ce ..+.C_ at each step in the division process to 


O 


the most recent values available for 


E ay 
S Da s m) 
dY | 
C =f (У...) 
n=l n-1l dt ( 38) 
_ ау 
Са = ТО o 


22 





In general 


Di 


С di ) 


ДЕ: (39) 


= г, ( 


о words, in the first step, the values of Ci are 


ЧУ 
шон 


result of the first division is a new value of Y to be used 


obtained from initial conditions for Y The 


in calculating Cn» EP T 2- for the second step in the 
division. New values to be used for the derivatives of Y 
with respect to t are calculated from the new values of the 
dependent variables. The process is difficult to describe 
MO ds DL ls simpre an concept. 


To illustrate the procedure an example is given: 


Example 3: 


ЧУ 2 
Let ЗЕ ү т (40) 
Carrying out the procedure outlined above, eq.(10) is 


written 


ЧУ 

— T е = 

dt Cay Ш (um 
The Laplace transform, assuming C to be constant and 


expressing in powers of 571 15 


Y(s) = LE (42) 
1 +Cs 


23 





Substituting the z-forms of Table I and dividing by T yields: 


ery a mare He) —^— 2 — (43) 
a (12 + 6ст) - 2441 + (12 -6ста2 
Letting T = 0.1, 
Y) - 0.1 + 1.0271 + 27° 
š 2 ee 
Carrying out the long division procedure yields: 
008 -1 -2 95273 zs 
КООЗ ТО 1002 +0. 1092 +#0.2052°2+0,3002 +... (4H) 


24 
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eu 





Е 322 033 126 223 010 
Figure 6 
К сае = 2,0 units/inch = Precise solution 
Y-scale = .5 units/inch + Approximate solution 
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TI l SUIS 


A. INTRODUCTION 
In order to test the present technique, a collection 
of linear differential equations with varying coefficients 


were studied, whose basic format is: 
Y + P(t)Y = Q(t) 


where P(t) and Q(t) are time functions. For the purpose of 
this work, nine equations in all were processed successfully 
for what P(t) took the values Vt, C and t^, and Q(t) took 

the values 1-u(t), + апа t^. 


The collection of equations, the case nubmers, and the 


sections in this work, are as follows: 


у +/ЕУ = 1 Сазе 1 
Yoa Бе =) Case 2 ТЕТЕ 
Y 4 t^Y 2 ] Case 3 
y «4t y = t Case ! 
Y T ҮҮ = t Case 5 J Tp UO 
Y + t^Y = t бэсе б 
ү + ЈЕ У = t? Case 7 
Y + tY = t2 Case 8 ION 
D t^y = Le Case 9 


In Case 10, an example from Ref. 2 was studied, Comments 


on this problem are given in Case 10 discussion. 


26 





In cases 1 to 9 the same type of procedure was followed, 
applying directly the steps as in II A of this work and 
experimenting for time intervals form 0.2 sec. to 2.0 sec. 
in order to show the advantages of this technique with 
relatively large At's, its stabllity together with an 
alsorithm for computer solution and application to hand 
calculation, 1.e., small number of division steps for good 


results. 


B. САЗЕ 1 — У +УМЕУ = 1 
The equation ae P(t)Y = Q(t) is rewritten for this 


саве ав: 


Sl 
KK 


ct 


+/ty = 1 1 (Еко) = 10 


| 
+ 
< 
1 
= 
hg 
n 
FI 


the Laplace transform is 


sy + Py = 1/s 


rearranging and expressing in powers of 571: 
s^y + Psy = 1 y (s^4Ps) = 1 
y 2 mi 


А 





Substituting the z-form values for ВК from Table I (page 15) 


Teme dose те 
2 DEDE 
A E 
Е 491 21-27. 
14 Р 2 
= 
1-7 


Dividing by T and rearranging: 


1 2 


mM+10Tz 4724 


n a (45) 
ОЕ Е Ср в 


Уд (2) = 


where again P is the time varying coefficient of Y in the 
differential equation and T is the iteration time, or in 
other words the time difference between results. 

The following step consists of dividing the denominator 
ме не штета ток апа substituting in P, its value for the 
corresponding time (t = nT). As the exponent of z for the 
first value in the quotient is zero, P(t) = ft at t = O, 
Pit = О in the first division step. 

The value for T first chosen was 0.2 sec. 

The algorithm used for computing the quotient (Program l, 
page 30) as the solution, for values of T from 0.2 sec. up 
to 2 sec. with the "precise" solution in continuous line and 
with the approximate solution as crosses follow; the "precise" 
solution was obtained by computer, using the subroutine 


called INTEG1. 
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Scales in the graphs (Figs. 7 to 20) are the same in 
раси сазе. 

In the numerical output tables, the time is referred to 
as T, the approximate solution is referred to as Q(N) and 
the "precise" solution as PR. TFIN is the final time. 

It can be seen that the approximated solution even for 


large time iterations is very close to the precise solution. 


FOR DELTA T = 0.20 TFIN = 5.2 
ELT Q(N) PR 
0.0 0): (8 ал (0 
0.2 ТӘЛІП 0.19300 
0.4 0.36010 0.36197 
076 0.49870 0.50039 
0.8 0.60588 0 060716 
20 0.682335 0.6841 
1.2 0.73455 0.73486 
je 0.76372 0575 36.3 
Re ӨЛШЕДІ? 02221211 
jS 0.77345 0.77289 
2.0 0 olo. 0 761205 
2,2 0.74384 0.74318 
2.1 0 70170 0.72109 
DS 0.69713 0.69689 
ао 02212215 0.67200 
3.0 0.64773 0.64738 
3.2 0.62392 0.62366 
3.8 ШЕСІ 0.60122 
30 0.58039 0.58027 
Bac 0.56092 0.56085 
lm 0.54297 0.54294 
4.2 0.52648 02 2011 
4.4 0.51132 0 "51153 
4.6 0.49738 От 
19: 0.48454 0.48458 
5.0 0.47268 0747272 
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1. Case 2 — Y + tY = 1 

The calculations are the same as in Case 1 with the 
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As in Case 1 and Case 2, the same calculations are 
followed and only in the division, the substitution of 


P(t) = tf 


is the sole difference. 
Output for several time iterations follows as in the 
BNO LS), 


previous cases (Figs. Computer algorithm is shown 


in Program 3 (page 63). 
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Equation (45) is repeated here for convenience: 
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C. CASE 1 – Ү +VtY =t 
For this case, in the equation Y+P(t)Y =.Q(t), 
P(t) =/t and Q(t) =t, so following the procedure the 


equation is rewritten as: 


ау 





cm cs = = 
ЗЕ РҮ t ша 0 
where P = /t . The Laplace transform is: 
= 2 
sy + Ру = 1/3 
so 
в Зу + Ps^y = 1 y (s? + Ps?) = 1 
53 + Рв“ 1 + Ps "Í 


Substituting the z-forms of 575 from Table I one gets: 


E 
Ya z) Ee аа a 
РЕЛЕ 
с 1-72 
Dividing by T and ordering in crescent powers of 277 
in numerator and denominator: 
Br. RE к. : 
Уд со = PEERS с зп суш 0500 з (46) 
(2+Pt) - (6+PT)zZ + (6-PT)z + (PT-2)z 
78 





Now, substituting for the time iteration wanted and 
replacing P = Vt for its value at the division step, the 
solution is obtained. 

In the same way as before, plots and outputs follow 
for time iterations from 0.2 sec. to 2.0 sec. (Figs. 49 to 
60). 

Algorithm for computer solution of this problem is shown 


in Program 4 (page 80). 
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ШЕ Саде 3 = ү + tY = t^ 


Same calculations as in case 7 with exception for 
the substitution on the value of P during the division 
steps, of Equation (47), repeated here, which now has the 
value P = t. 

The algorithm for computer solution (Program 8, 
page 144), the output for several time iterations and the 


respective graphs follow (Figs. 104 to 117) 
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8. Case 9 — Y + t Y = t^ 


Same reasoning as in case 8. This time P = t°. 
Computer algorithm for solution of this problem ís 
shown in Program 9 (Page 160). Solutions and graphs for 


Several T's follow in Figures 118 to 131. 
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Equation (47) is repeated for convenience: 
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E. CASE 10 
No one of the previous cases has a time varying 
coefficient which is exponential. So, from Ref. 2 this 


case was extracted. The equation is: 
de e С 
Б = = nn 
С dt a 0 being r(t) R exp (7 


where Т. is a factor in the variation of r. Dividing all 


the terms by C yields: 


de 1 -t 
p EC CN) 
dt То : Т. 

where T5 = ВС. 


For this example T, and T, were made equal to one and 
jJ 2 


= 1.0 . Rewriting the equation and expressing its 
1; 


“(t=0) 
Laplace transform in powers of s. 


e + Ce = O C = e 


5 РСЕ Е 0 


1 1 Св. 


псе rules in Example 2, Section ІЛЕ were followed, 


k 


substituting the z-forms of s ”, dividing by T and adding 


enema ccults of the Quotient to 1,0 интер is the time 


function for l/s: 


175 





Ще ре rer 
па A "un TE 
Е Т 1421 
l-z 
bor T - 0,1 
T О O 


CP A A O a 


Dividing the denominator into the numerator and adding one 


to the result one gets 
0.99206 + 0.91412571 + 0. 8500777“ + 
However, looking at the graph (Fig. 132) the precise 


and the approximated solution start diverging after a while. 


mars 1S due to the fact that in 





ис сое vo zero the equation becomes 


E = 1/5 


шан is the final value that the approximated solution 


approaches. 
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This fact is very important because it leads to these 
important conclusions: 

T. If Cy goes to zero with time, first check the Laplace 
transform and the initial equation substituting Cy by zero. 
If by inspection they do not lead to the same type of 
solution, this method does not apply. 

2. Verify, before doing the division if any of the 
varying coefficients, unexpectedly, goes to infinity; by 


Gividing by zero, for instance. 
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IV. CONCLUSIONS 


This technique, for the successfully area covered in 
this paper (linear differential equations with time varying 
coefficients), showed to have a good precision even for 
relatively large iteration steps. The labor necessary to 
complete a problem by handwork is not very involved and 
needs only a small algorithm for computer solution with 
C.P.U. times much smaller than existing numerical methods. 

However, with time varying coefficients that go to zero 
with time, attention must be paid to the application of 
this method. As a general rule, this method does not apply. 

In some cases, the final equation for the approximate 
solution, before the divisicn is performed, must be altered 
in order to not have time varying coefficients (referred as 
C, in Section IIc) in the numerator. Caution is necessary 
to verify, before starting the division steps if any varying 
coefficient Ва has a value of infinity; by division 
by zero when substituting a particular value of Cy > 

The technique tested in this paper does not apply to the 
case of linear differential equations with varying coefficients 
only. Linear differential equations with constant coefficients 
and nonlinear differential equations are, as well, solved by 
this technique. 

These areas, in the opinion of the author, are valuable 


as Coples lOr other theses. 
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